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ABSTRACT 


A large set of numerical simulations of magnetohydrodynamic (MHD) turbu¬ 
lence induced by the magnetorotation al instability ('MRI ') is presented. Revisit¬ 
ing the previous survey conducted bv ISano et al.l (120041) . we investigate the gas 
pressure dependence of the saturation level. In ideal MHD simulations, the gas 
pressure dependence is found to be very sensitive to the choice of a numerical 
scheme. This is because the numerical magnetic Prandtl number varies accord¬ 
ing to the scheme as well as the pressure, which considerably affects the results. 
The saturation level is more sensitive to the numerical magnetic Prandtl number 
than the pressure. In MHD simulations with explicit viscosity and resistivity, 
the saturation level increases with the physical magnetic Prandtl number, and 
it is almost independent of the gas pressure when the magnetic Prandtl number 
is constant. This is indicative of the incompressible turbulence saturated by the 
secondary tearing instability. 


Subject headings: instabilities — magnetohydrodynamics (MHD) — turbulence 
— methods: numerical 


INTRODUCTION 


Space is hlled with dilute, magnetized plasmas. Plasmas that exhibit sufficiently large 
Reynolds and magnetic Reynolds numbers are expected to be in a turbulent state. For 
example, weakly magnetized accretion dis ks are subject to magiietic tu rbulence induced 
by the magnetorotational instability (MRI; iBalbus Sz Hawlevlll99ll. Il998l ). Since numerical 
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simulations have revealed that the Maxwell stress in the MRI-induced turbulence is orders 
of magnitude larger than the viscous stress estimated from ordinary molecular viscosity, the 
MRI-induced turbulence is believed to play a critical role in the o utward angular momentum 
transport and subsequent mass accretion toward the central star flHawlev et ahlll995l. Il996l) . 
The nonlinear stage of turbulence in magnetized plasmas has been extensively studied via 
magnetohydrodynamic (MHD) simulations. 

Ideal MHD equations without physical viscosity and resistivity are often utilized for 
the numerical study of astrophysical plasma dynamics. Artificial viscosity and resistivity 
sh ould be explicitly added t o central (e.g., Lax-Wendroff) or operator-splitting (e.g., ZEUS 
by Stone fc Norman f 1992alf3) I schemes to suppress spurious oscillations and capture shocks. 
On the other hand, Godunov-type schemes automatically embed sufficient numerical diffu¬ 
sion because of their upwind property by virtue of a Rieman n solver. Recent God unov-type 
schemes such as NIRV ANA f Ziegler 2004 . 2008), RAMSE S ( Fromang et ah 2006 ). PLUTO 
( Mignone et ah 2007 ). and ATHENA ( Stone et al. 2008 ) have been designed to capture 
shocks and discontinuities via exact or approximate Riemann solvers, and to resolve small- 
scale turbulence via high-order reconstruction techniques. The spatial accuracies of the 
fluid and magnetic held in Godunov-type schemes are not necessarily identical because of 
the treatment of the induction equation (divergence or curl form) and spe c ial co nsideration 
given to the divergence-free condition for the magnetic held. iKritsuk et al.l (1201 if) have com¬ 
pared the accuracy of nine MHD simulation codes from a turbulence decay simulation, and 
they have shown that a code with high accuracy for velocity does not necessarily possess 
high accuracy for magnetic held, and vice versa. The numerical magnetic Prandtl number 
that is dehned as the ratio of the numerical viscosity to resistivity varies according to the 
applied schemes, and it can be smaller or larger than unity. 

The magnetic Prandtl number Pr^ determines and controls essential properties of MHD 
turbulence. The num ber relates the viscous dis sipation scale ly to the resistive dissipation 


scale In- For example, ISchekochihin et al.l (120041) have argued that the kinetic energy input 


at the system scale cascades down to the Kolmogorov dissipation scale /j,, and subsequently, 
viscous eddies amplify small-scale magnetic helds down to ~ Prm ' G for Pr^ > 1. When 
Spitzer values are used for the viscosity and resistivity of fully ionized, weakly magnetized 
plasmas, the magnetic Prandtl number is estimated to be Pr^ 10 where T an d 

n denote the temperature and number density in cgs units, respectively (jSpitzerl 119621) . 
Astrophysical objects can have magnetic Prandtl numbers that lie far from Pr^i = 1. There 
are extreme environments that have Pr^ S> 1 (e.g., stellar corona, active galactic nuclei 
disks) or Pr^ -C 1 (e.g., stellar interior, young stellar object disks). In such cases, there is 
a very wide sub-viscous or sub-resistive range over which the fluid and magnetic held are 
decoupled with each other. It can considerably alter the generation and dissipation processes 
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of MHD turbulence with respect to those observed for Prm = 1. ISchekochihin et ahl (120041) 
have performed numerical simulations of forced MHD turbulence, and they have shown that 
a small-scale dynamo with S> 1 generates magnetic helds in the sub-viscous range. The 
magnetic held energy positively correlates with Pr^ and it does not satisfy scale-by-scale 
equipartition with the kinetic energy in the saturated state. Such physics is not included in 
ideal MHD simulations, whose numerical magnetic Prandtl number is of the order of unity. 

When ideal MHD equations are solved with a Godunov-type scheme to study MHD 
turbulence, the viscous and resistive dissipation scales and hence the magnetic Prandtl num¬ 
ber are implicitly determined depending on the choice of the reconstruction technique and 
the Riemann solver. It is of critical importance to clarify effects of numerical dissipation 
on practical MHD turbulence problems. To this end, we conduct a large set of ideal MHD 
simulations of the MRl. We revisit an extensive survey of the gas pressure depende nce of 
the saturation level of the MRl-induced turbulence as conducted bv ISano et ahl (120041) . Our 
simulation reveals that the saturation level and its pressure dependence are very sensitive to 
the numerical magnetic Prandtl number, thereby indicating the necessity of explicit viscosity 
and resistivity to control the magnetic Prandtl number. Consequently, we conduct visco- 
resistive MHD simulations in order to investigate the parameter dependence and discuss the 
saturation mechanism of the MRl-induced turbulence. 


This paper is organized as follows. In Section |2l we describe our numerical method based 
on a Godunov-type scheme. In Section El we present simulation results. First, we assess 
the effect of numerical dissipation on the saturation level of the MRl-induced turbulence 
with the ideal MHD simulation in Section 13.11 Subsequently, we investigate the parameter 
dependence of the saturation level with the visco-resistive MHD simulation in Section 13.21 
Section 0] discusses the saturation mechanism of the MRl-induced turbulence based on the 
simulation results. Finally, we summarize the paper in Sectional 


2. NUMERICAL METHOD 


We perform th ree-dimensional cornpressible MHD simulati ons under a local shearing 
box approximation (iHawlev et al.lll995l: IStone fc Gardiner! 120101) . Without considering the 
viscosity, resistivity, and gravity stratihcation, the relevant equations are as follows. 


|^ + V.(pu)=0, 

dpu' 


dt 


puu 


^ BB 

P + — 1- 

Svr / dvr 


= —2H X pu' + gUpu^Ey, 


( 1 ) 

( 2 ) 
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( 3 ) 

( 4 ) 


dB 

dt 

E = 


+ W X E = 
—u X B, 


0 , 


where p denotes the density, P the gas pressure, B and E the magnetic and electric helds, 
u = uq + u' the total velocity, Uq = —qQxBy the background flow in the shearing box 
(local linear expansion of the rotating motion), u' the deviation from the background shear 
flow, fl = Qbz the rotation axis, q the shear parameter {q = 3/2 for Keplerian disks is 
used). The entropy equation is omitted because we adopt an isothermal equation of state. 
P/p = cl = constant where Cg is the sound speed. 


2.1. Numerical Algorithm 


The equations are advanced based on a technique proposed bv IStone fc Gardinerl (120101) . 
Their technique decomposes the equations into the MHD part including the shearing source 
terms (right-hand-side of equation ([2])) and the orbital advection part with the background 
shear flow Uq. The latter part is advanced as simple one-dimensional linear advection equa¬ 
tions because uq has only a |/-component and is independent of y and time. The shearing 
source terms are integrated by means of a second-order Crank—Nicholson method so as to 
guarantee the conservation of epicyclic motion in a discrete sense. 

For the MHD part, we adopt a Godunov-type scheme (Minoshima, T., Matsumoto, 
Y., and Miyoshi, T., in preparation). The scheme involves (i) nonlinear reconstruction of 
characteristic variables from cell centers to cell faces, (ii) evaluation of numerical fluxes 
dimension-by-dimension with a Riemann solver, and (hi) update of the induction equation 
with a method similar to the upwind constrained transpor t method flLondrillo fc Del Zanna 


2000l: iLondrillo fc del Zannal |2004 iDel Zanna et ahl 120071) . The induction equation is ad¬ 


vanced in curl form to preserve the divergence-free condition for the magnetic held in a 
discrete form to machine accuracy. Details of this method are described in Appendix 
For nonlinear reconstruction, we adopt a fifth-order Weighted-EN O scheme (WENO ; 


Jiang &: Shulll996l) or a combination of a fifth-order WENO-Z scheme flBorges et al.ll2008l) 


and a monotonicity-preserving scheme (jSuresh fc Hnvnhl 119971: iBalsara fc Shul 120001) . Here¬ 
after, the latter is termed as the WZMP scheme, whose accuracy is slightly better than that 
of the WENO scheme. For Riemann sol vers, we adopt the single-state Harten—Lax—van 
Leer (HLL) approximate Ri emann solver (IHarten et al.l 119831) or the multi-state HLLR ap¬ 


proximate Riemann solver ( Mivoshi fc Kusa'ii^ 


200 


The former solver resolves the fast 


mode, and the latter resolves the Alfven mode as well. 


As is described in Appendix our simulation code can use different reconstruction 
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functions and Riemann solvers for the fluid and magnetic field. Consequently, we refer to 
the applied scheme as F1-R1-F2-R2, where “FI” and “Rl” denote the reconstruction func¬ 
tion and Riemann solver for the fluid, and “F2” and “R2” denote those for the magnetic held, 
respectively. For example, the WENO-HLLR-WZMP-HLL scheme uses the WENO recon¬ 
struction and the HLLR Riemann solver for the huid update, and the WZMP reconstruction 
and the HLL Riemann solver for the magnetic held update. By varying the combination 
of the reconstruction function and the Riemann solver for the huid and magnetic held, we 
assess the ehect of numerical dissipation on the MRI-induced turbulence. 


We adopt a hnite volume approach. Finite volume schemes require reconstruction of the 
numerical hux over c ell faces to achieve a desired order of accuracy in multidimension (e.g., 
Titarev fc Toroll2004l) . We ignore this reconstruction for simplicity, that is, we approximate 
F = F{U) where U, F are the physical variable and the corresponding hux, and the overline 
represents the average over the cell face. This approximation is inaccurate when E' is a 
nonlinear function of U. It degrades the order of accuracy in multidimension even using a 
high order reconstruction dimension-by-dimension. 

We examine the order of accuracy of our scheme from the propagation of circular- 
polarized Alfven waves in two-dimensional homogeneous medium. The propagation angle is 
45° with respect to the x axis and the plasma beta value is 0.1. We use the WENO-HLLR- 
WEN O-HLL scheme and the third-order TVD Runge-Kutta time integration fjShu Sz Osher 
19881) . The order of accuracy of the scheme is presented in Figured] The scheme preserves 
hfth order of accuracy at a coarse grid size, but the accuracy becomes worse than the fourth 
order at a hne grid size. In this problem, the order of accuracy can be retrieved by increasing 
the order of the Runge-Kutta time integration, implying that multi-dimensional properties 
are somewhat incorporated by using the multi-step time integration. The similar result is 
obtained with the WZMP-HLLR-WZMP-HLL scheme. However, this does not necessarily 
hold for practical problems in highly inhomogeneous medium. The order of accuracy may 
be degraded as much as the second order. In following simulations, we use the third-order 
TVD Runge-Kutta time integration. The fourth-order time integration does not make a 
signihcant difference. 


Stone fc Gardiner! (120101 ) have pointed out that the same high-order reconstruction func¬ 


tion used in the computational domain is required for the remapping of variables to the 
shearing boundary to maintain the accuracy of the calculation. We conhrm that a simple 
arithmetic average at the shearing boundary severely degrades the accuracy of the simu¬ 
lation. Consequently, we remap the variables by solving the linear advection as follows. 
The shearing distance at the left boundary relative to the right, Lg = qQLt, is divided into 
integer and non-integer numbers of grids, N = int(Ls/A) and D = Lg/A — N. Subse- 
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quently, the variables at the right boundary are remapped to the left boundary such 
that U^{yj) = U^{yj-N — DA). We use the same reconstruction function to shift with 
a non-integer distance of DA. 


The growth of the MRI is very sensitive to the amplitude of the net magnetic flux. 
The use of the nonlinear rec onstruction and the sh earing periodic boundary violates the 
conservation of the variables flGressel fc Ziegler! 120071) . To avoid the problem, we correct the 
numerical flux of (that is, E.^,) at the she aring boundary (x = ±0.5) in a manner similar 
to that adopted by IStone fc Gardinerl (120lOl) . 


2.2. Initial Setup 


The space, time, and velocity are normalized by the system size L, the rotation period 
and the rotation velocity LQ, respectively. The simulation domain of [—0.5, 0.5] x 
[—2,2] X [—0.5, 0.5] is resolved with grid points of 32 x 128 x 32 (grid width A = 1/32) in 
Section [TT] High-resolution simulations with grid points of 64 x 256 x 64 (A = 1/64) are 
also conducted to check numerical convergence in Section [3A2] and include explicit viscosity 
and resistivity in Section 13.21 The boundary condition is periodic along the y (azimuthal) 
and 2 ; (vertical) directions, and shearing periodic along the x (radial) direction, respectively. 
The initial condition is uniform, i.e., {p,u'^,u'y,u'^, P, B^, By, 3^) = (1, 0, 0, 0, Pqj 0, 0, Hq), 
and we impose an incompressible small random noise to uniformly in the y — z plane 
and constantly in the x direction to initiate the instability. The parameter Bq is hxed to be 
0.025 throughout the study, and the corresponding wavelength of the fastest growing mode 
(FGM) of the MRI is A t^cm = = 27rn^/f2 = 0.16 where va = Bq/^ denotes the Alfven velocity. 
Following ISano et al.l (120041) . we investigate the gas pressure dependence of the saturation 
level by varying the initial plasma beta value as 10^, 10^, 10^, 10^, and 10®. The simulation is 
relevant to actual accretion disks when the vertical scale size L is comparable to the pressure 
scale height \/‘lPjp/D at /3 = 10^. The simulation runs more than a 100 orbit periods to 
improve the statistics. We measure the statistical variables by averaging in whole space and 
time over the last 50 orbit periods. 


3. RESULTS 

3.1. Ideal Magnetohydrodynamic Simulations 

Tables [U to [3] list the simulation parameters and results obtained for the ideal MHD 
simulation. As can be noted from the tables, we use different Riemann solvers for the fluid 
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and magnetic field updates. Table [T] lists results for the application of the HLLR Riemann 
solver for the fluid and the HLL Riemann solver for the magnetic held. Table |2] lists the 
results for the case when the HLL Riemann solver is used for both the huid and magnetic 
held. Table [3] corresponds to the case when the HLLR Riemann solver is used for both 
the huid and magnetic held. Column (1) specihes the model applied. The hrst letter I 
denotes “ideal”. The second letter denotes the logarithmic of the plasma beta value (listed 
in column (4)). The remaining letters represent abbreviation of the reconstruction function 
and the Riemann solver (W = WENO, Z = WZMP, H = HLL, R = HLLR). Columns (2) 
and (3) specify the reconstruction function used for the huid and magnetic held, respectively. 
Columns (4) and (5) list the initial plasma beta and gas pressure, respectively. Columns 
(6) and (7) list the statistical average of the magnetic energy Em = /Sir and the kinetic 
energy = pu'‘^/2, respectively. Columns (8) and (9) list the statistical average of the 
Maxwell stress wm = a nd the Reynolds stress wr = pu'^Uy, respectively. Column 

(10) lists the a parameter given bv IShakura &: SunvaevI (119731) . a = {wm + wr)/Pq. Finally, 
column (11) lists the “numerical” magnetic Prandtl number, which is estimated in a later 
section. 

First, we verify the correlation of statistical averages in the turbulent state (Figure 
[2]). All the simulations listed in Tables [I]l3] are considered. The Maxwell stress is smaller 
than the magnetic energy by a factor of ~ 2 and larger than the Reynolds stress by a 
factor of ~ 5 within the range of two orders of magnitude. The Maxwell stress, magnetic 
energy, and the Reynolds stress sh ow good correlation , and the results appea r consistent 


with those of previous studies (e.g., iHawlev et aLlll995L Il996l: ISano et al.l 120041) . Figure [3] 
shows time prohles of the Maxwell and Reynolds stresses in three runs denoted as (a) 12- 
WRZH (Po = 0.03125), (b) I4-WRZH (Pq = 3.125), and (c) I6-WRZH (Pq = 312.5). Many 
random spikes are prominent for both the Maxwell and Reynolds stresses in the turbulent 
state for the runs I4-WRZH and I6-WRZH. Their amplitude is comparable to or larger than 
the temporally-averaged values. When the pressure is low (I2-WRZH), on the other hand, 
the temporally-averaged value is smaller than that in the high pressure cases and spikes are 
less prominent. Figure H] shows the ^/-component of the magnetic held at the (a) peak and 
(b) decay of a spike in the run I4-WRZH. Large-scale hows are well developed at the peak. 
The amplitude of the hows is comparable to the background shear how. The magnetic held 
is stretched and amplihed predominantly along the ^/-direction. The spatially-averaged 
magnetic energy is increased by a factor of several hundred from the initial value, thereby 
giving plasma beta of several tens. Consequently, the how velocity that is comparable to the 
Alfven speed is still much slower than the sound speed. At the decay, the ho w collapses to 
a smaller scale due to secondary parasitic instabilities flGoodman fc Xul 119941) . The growth 
of the large-scale how and its subsequent decay repeatedly appear in the turbulent state. 

















These features are in good agreement with the results obtained bv ISano fc Inutsukal (120011) . 
Therefore, we can conclude that the simulation successfully solves the MHD turbulence 
induced by the MRI. 


Following ISano et al.l (1200411 . we investigate the gas pressure dependence of the satura¬ 
tion level of the Maxwell stress in Figure [5])a). The data are taken from Tabled] that is, we 
use the HLLR Riemann solver for the fluid and the HLL Riemann solver for the magnetic 
held. There are three combinations of the reconstruction function for the huid and magnetic 
held. The crosses correspond to the results obtained with the WENO reconstruction for both 
the huid and magnetic held (I?-WRWH)0. The asterisks represent the results of the WENO 
reconstruction for the huid and the WZMP reconstruction for the magnetic held (I7-WRZH). 
This model improves the accuracy of the magnetic held with respect to the hrst I7-WRWH 
model. The diamonds denote the results obtained with the WZMP reconstruction for both 
the huid and magnetic held (I7-ZRZH). This model improves the accuracy of the huid as 
well as the magnetic held with respect to the hrst I7-WRWH model. 


Sano et al.l (120041) have shown that the saturation level of the Maxwell stress weakly 
depends on the gas pressure as wm oc for the isothermal case. However, the pressure 
dependence in our simulation is considerably diherent from their result. In the low-pressure 
range ofP = 0.01 — 1.0, the Maxwell stress increases with increasing the pressure regardless 
of the choice of models. With increase in the range of pressure P = 1.0 — 1000, the pressure 
dependence dihers among the models. In the I7-WRWH model, the Maxwell stress is a 
decreasing function of the pressure. In the I7-WRZH model, the Maxwell stress is a weakly 
increasing function of the pressure. In the I7-ZRZH model, the Maxwell stress is indepen¬ 
dent of or weakly decreases with the pressure. The difference is more signihcant at higher 
pressures. Figure [5])b) shows the time prohle of the Maxwell stress at the highest pressure 
with the two models I6-WRWH and I6-WRZH. There are signihcant differences in not only 
the average but also the peak and huctuation amplitude between the models. This result 
indicates that our simulation is very sensitive to the reconstruction function, and hence, it 
fails to evaluate the gas pressure dependence. 


We also examine the effect of the difference of the the Riemann solver on the gas pressure 
dependence of the saturation level. In the previous case, there is a mismatch of the Riemann 
solver between the huid (HLLR scheme) and the magnetic held (HLL scheme). Here, we 
use the same Riemann solver for both the huid and magnetic held. We hrst use the HLL 
scheme, which degrades the accuracy of the huid with respect to that in the previous case. 
Secondly, we use the HLLR scheme, which improves the accuracy of the magnetic held from 


^The symbol “?” expresses an arbitrary digit from 2 to 6. 
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the previous case. The results are listed in Tables |2] and [3l and the pressure dependence 
of the saturation level of the Maxwell stress is shown in Figure [6l These results appear to 
be less sensitive to the reconstruction function than in the previous case, thereby indicating 
that the use of the same Riemann solver for both the fluid and magnetic held improves 
the convergence of the simulation. However, the results are found to be very sensitive to 
the Riemann solver itself. The Maxwell stress is almost independent of the pressure in the 
high-pressure range of P = 1.0 — 1000 with the F1-HLL-F2-HLL scheme (Figure [6]^a)). The 
amplitude is comparable to or slightly larger than that obtained with the F1-HLLR-F2-HLL 
scheme. On the other hand, the Maxwell stress is an increasing function of the pressure 
within the explored parameter range with the F1-HLLR-F2-HLLR scheme (Figure in](b)). 
The amplitude is considerably larger than those in the previous cases with the Fl-HLLR- 
F2-HLL and F1-HLL-F2-HLL schemes, particularly at higher pressures. I n the intermediate 
pressure range of P = 0.1 — 100, the result follows the correlation by ISano et al.l (120041) 
(dash-dotted line). 


The effects of numerical dissipation on the gas pressure dependence of the saturation 
level are briefly summarized as follows: (i) When the accuracy of the magnetic held is 
improved, the saturation level increases, (ii) When the accuracy of the huid is improved, 
the saturation level decreases. This trend indicates that the saturation level depends on the 
numerical magnetic Prandtl number, Pr^ = where u,r] are the kinematic viscosity and 
resistivity, respectively. 


3.1.1. Numerical Magnetic Prandtl Number 


The m agnetic Prandtl number d ependence of the MRI-induced turbulence was hrst 
reported by Lesur fc Longaretti J2007 ) for fi nite net magnetic hux cases (P^ = const), and 
it has been studied bv iFromang et al.l (120071) for zero net hux cases oc sin(27ra;)). They 
have carried out three-dimensional visco-resistive MHD simulations under the shearing-box 
approximation, and they have shown that the satur ation level of the Maxwel l stress is an 
increasing function of the magnetic Prandtl number. Lesur fc Longarettil (200^^ have shown 


that wm oc for the explored range of 0.12 < Pr^ < 8.0. IFromang et al.l (120071) 

have determined the critical magnetic Prandtl number as Pr^,ait ~ 1 for zero net hux cases, 
below which turbulence is quenched. 

Our ideal MHD simulations could not evaluate the gas pressure dependence of the 
saturation level of the M axwell stress. It depend s on the choice of the num erical scheme. 


Based on the studies bv iLesur fc Longarettil (l2007l) and IFromang et al.l (120071) . we speculate 
that the numerical magnetic Prandtl number depends on the choice of the numerical scheme 
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and varies with the pressure. To verify it, let us consider the propagation of incompressible 
Alfven waves along the x-direction, 


dvt 

dBt 

dt 

4:71 p dx 

dBt 

_ p dvt 

dt 

dx 


(5) 

( 6 ) 


where Vt and Bt denote the transverse components of the velocity and magnetic held, p and 
> 0 are uniform, and = 0 is assumed. As is well known, the combination of these 
equations gives the advection equation for the Elsasser variables ± VAt, 


9/ A 

Oj. T "^Ax o — 0. 

at ox 


(7) 


To evaluate the numerical error in the Alfven wave propagation problem, we discretize 
the equations ([5]) and ([6]) in space by means of a hrst-order local Lax—Friedrichs scheme 
(equivalent to the HLL scheme in the homogeneous case). 


dvt,i 

dt 

dBt,i 

dt 


'^Ax 


\/47rp 


2A 




D I ^ z' o O Z3 I D 

-Ox- ^ -h ^ [Bt,i+i - ZBt^i + Bt,i-i 


i-l) + O (A) , 
) + O (A) , 


( 8 ) 

(9) 


where c denotes the maximum speed among eigenmodes. These equations yield the advection 
equation for the Elsasser variables in a discrete form as 


dt 

df~ 


'^Ax 


+ VAx 


ft +1 - ft ^ \c-VAx\Aft+i-2ft + fU 
A 2 A2 

/h - /Ai _ \c-VAx\^ /i +1 - 2/h + /Ai 


+ 0(A), 
+ 0(A). 


( 10 ) 

( 11 ) 


dt A 2 A2 

The left-hand-side corresponds to the hrst-order upwind discretization. In addition to the 
discretization error 0(A), there is a second-order dihusion term on the right-hand-side, which 
acts as the resistivity to dissipate Alfven waves. The coefficient of this numerical dihusion, 
^num = |c — vax\^/2, is negligible for low beta (c = Cfast ~ vax), but not for high beta, 
^num oc Cfast oc The use of single-state Riemann solvers such as the Lax—Friedrichs and 

HLL schemes leads to the pressure-depe ndent num erical dihusion of Alf ven waves. Multi- 
state Riemann solvers such as the Roe fjRoelll98lh . HLLR, and HLLD flMiyoshi fc Knsano 


20051 ) schemes signihcantly reduce the amount of numerical dissipation because they are 


designed to resolve the Alfven mode as well as the fast mode. 


However, the above statement holds only for one-dimensional simulations. In actual 
multi-dimensional simulations, even multi-state Riemann solvers lead to the pressure-dependent 
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numerical diffusion of Alfven waves owing to an error associated with the directional split¬ 
ting (the one-dimensional Riemann solver is applied dimension-by-dimension). Numerical 
diffusion of MHD turbulence is inevitable, and it depends on the plasma beta. Therefore, we 
consider that the pressure dependence of the saturation level of the MRI-induced turbulence 
is indistinguishable from the numerical magnetic Prandtl number dependence. 


To quantitatively discuss the effect of the numerical magnetic Prandtl number on the 
saturation level of the Maxwell stress, we estimate numerical viscosity and resistivity as fol¬ 
lows. We restart simulations with = 0 to numerically dissipate turbulence. Subsequently, 
the spatially averaged decay rates are measured as the numerical viscosity i^num and resis¬ 
tivity r^num via the following equations obtained by integrating viscous momentum equation 
and resistive induction equation over space. 


^num 


Vrmm 


d n'2 
dt 2 


f W ■ {j X B)\ 



/ 

(V X u'f + (4/3) (V • u'f 

1 P ij 




d 

dt Svr 


■ {j X B) 



( 12 ) 

(13) 


where j = V x B/dvr denotes the current density, and the overline represents spatial aver¬ 
aging. The hrst term in the numerator represents the energy change, and the second term 
represents the work done by the dynamo. The spatial variation in the density is ignored for 
simplicity. Note that this method is inaccurate for strongly compressible cases. Finally, we 
estimate the numerical magnetic Prandtl number as Prm,num = ^num/hnum- This is a crude 
estimation. In reality, numerical diffusion coefficients depend on th e wavenumber, thus de¬ 
pend on the period over which the MHD simulation data are taken. iFromang fc Papaloizou 
(120071 ) have proposed a more rigorous treatment to estimate the numerical resistivity by an¬ 
alyzing the induction equation in Fourier space. We adopt the above method for simplicity, 
and we assume that the numerical viscosity and resistivity exhibit a similar dependence on 
the wavenumber, and consequently, the numerical magnetic Prandtl number is insensitive to 
the wavenumber. 


Example estimations of the numerical magnetic Prandtl number are shown in Figure 
[3 for the WENO-HLLR-WENO-HLL scheme. In this scheme, the HLL Riemann solver is 
used only for the magnetic held, which degrades the accuracy of the magnetic held when 
compared with that of the huid as the pressure increases. Therefore, the numerical mag¬ 
netic Prandtl number is expected to be a decreasing function of the pressure. This is con- 
hrmed in our estimation; the temporally averaged values are 0.76,0.43,0.51,0.45,0.36 for 
p = 102, lO^, 10^ 10^ 10®, respectively. 

The averaged values of the numerical magnetic Prandtl number in all simulation runs 
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are listed in column (11) of Tables [I]l3l We observe the following systematic trends: (i) 
When both the fluid and magnetic held are updated with the HLL scheme (Table [2]), the 
numerical magnetic Prandtl number is close to unity and almost independent of the pressure, 
(ii) When the accuracy of the huid is improved by the HLLR scheme (from Table |2]to Table 
[1]), the numerical magnetic Prandtl number clearly decreases and is a decreasing function 
of the pressure, (iii) When the accuracy of the magnetic held is improved by the HLLR 
scheme (from Table [T] to Table [3]), the numerical magnetic Prandtl number clearly increases, 
(iv) When both the huid and magnetic held are updated with the HLLR scheme (Table 
ED, the numerical magnetic Prandtl number is larger than unity, and it is an increasing 
function of the pressure, (v) When the WZMP reconstruction function is utilized to improve 
the accuracy of the magnetic held, the numerical magnetic Prandtl number clearly increases, 
(vi) When the WZMP reconstruction function is utilized to improve the accuracy of the huid, 
the numerical magnetic Prandtl number clearly decreases. The numerical magnetic Prandtl 
number is smaller or larger than unity, and its pressure dependence is not straightforward, 
instead depending on the choice of the Riemann solver as well as the reconstruction function. 

In Figure [H](a), the Maxwell stress is shown as a function of the numerical magnetic 
Prandtl number. Error bars show the standard deviation of the temporal average. Positive 
correlation is clearly seen between the two parameters. If we exclude the data in strongly 
compressible cases (/3 = 10^, triangles) and the data with extremely large errors, we obtain 
the correlation wm oc (dash-dotted line). In Figure IHl(b), the Maxwell stress is shown 

as a function of the gas pressure and the numerical magnetic Prandtl number. Except for 
strongly compressible cases, the Maxwell stress is more sensitive to the numerical magnetic 
Prandtl number than the pressure. Solutions of the ideal MHD simulation of the MRI- 
induced turbulence are found to be subject to the numerical magnetic Prandtl number of the 
applied scheme. Therefore, we conclude that it is necessary for the MRI-induced turbulence 
simulation to use explicit viscosity and resistivity to control the magnetic Prandtl number. 
The visco-resistive MHD simulation is indispensable to discuss the property of the MRI- 
induced turbulence. 


3.1.2. Convergence 

In ideal MHD simulations, numerical dissipation scales are related to numerical resolu¬ 
tion (number of grid points) as well as characteristics of the applied scheme. Previous studies 
have addressed the numerical convergence of solutions of the MRI under an unstratihed local 
shearing box approximation, and they have found that it depends on the initial magnetic 
held conhguration; increasing the resolution decreases the saturation level for zero net hux 
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cases flFromang et al.l 120071: IPessah et al. 
level for finite net flux caseJif Bodo et al. 


2007 


2008 


Silversll2008[) . 


Bodo et al.l 201 ih whereas it increases the 


To address the convergence of our simulation results, we conduct the same runs with 
doubling the resolution (grid points of 64 x 256 x 64). We use the three combinations of 
the Riemann solver that are same as the previous low resolution runs, and the WENO 
reconstruction is adopted for both the fluid and magnetic held throughout this subsection. 
Table 0] lists simulation parameters and results. The format is same as that used in Tables 
[T]12] except that the hrst letter H in column (1) denotes “high resolution” and columns (2) 
and (3) specify the Riemann solver used for the huid and magnetic held, respectively. 


Figure [9] compares the Maxwell stress in the low and high resolution runs. Except for 
the case at P = 312.5 with the WENO-HLLR-WENO-HLLR scheme, the Max well stresses 
i n the high resolution runs are larger than those in the low resolution runs. iGuan et al. 


(12009[) have argued that the increase in the saturation level is attributed to the resolution 
of small-scale structures near the correlation length at which energy is injected. We do not 
achieve the convergence within currently available computational resources. Meanwhile, we 
hnd that the numerical magnetic Prandtl number dependence are still observed in the high 
resolution runs; the Maxwell stress is low in the H7-WRWH model, medium in the H?- 
WHWH model, and high in the H7-WRWR model. Using the method in Section 13.1.11 the 
numerical magnetic Prandtl numbers are estimated to be Pr^,num ~ 0.5,1.0, and 2.0 in the 
H7-WRWH, H7-WHWH, and H7-WRWR models, respectively. The Maxwell stress roughly 
obeys the expression wm oc which is similar to that obtained in the low resolution 

runs. Note that the numerical magnetic Prandtl number is the ratio of the numerical viscosity 
to resistivity which are inversely proportional to the resolution of fluid and magnetic held, 
respectively. If the numerical scheme is designed so that the orders of accuracy of huid and 
magnetic held are nearly equal, the numerical magnetic Prandtl number is expected to be 
insensitive to the number of grid points. It rather sensitive to characteristics of the applied 
scheme. This implies that diherent ideal MHD simulation codes may lead to diherent results 
depending on the numerical magnetic Prandtl number, even if sufficient grid points are used 
to conhrm numerical convergence. 


^This will hold at low and moderate resolutions. Since it cannot grow indefinitely, the level should be 
converged to an asymptotic value for high enough resolution. 































14 


3.2. Visco-resistive Magnetohydrodynamic Simnlations 


We perform visco-resistive MHD simulations of the MRI to examine the gas pressure 
dependence of the saturation level at a constant magnetic Prandtl number. The viscous and 
resistive terms are advanced by an operator-splitting method to primitive variables, 


du' 

dt 

dB 

dt 


-V X V X n' -I- -V (V ■ u') 
3 


—rjV xV X B, 


(14) 

(15) 


where uniform viscosity and resistivity are assumed. A simple second-order central dif¬ 
ference method is used to discretize these equations. The viscosity and resistivity vary 
as (0.5,1,2,4) X 10“^, and the “physical” magnetic Prandtl number Pr^ = varies as 
0.5,1,2,4. The initial plasma beta value varies as 10^,10^,10^,10® (omitting the strongly 
compressible case oi(3 = 10^). The WENO reconstruction and the HLLR Riemann solver are 
used for both the fluid and magnetic held throughout the visco-resistive MHD simulation. 
The number of grid points is 64 x 256 x 64. The other parameters and numerical techniques 
are identical to those used in the previous ideal MHD simulation. 

Table [5] lists simulation parameters and resnlts obtained for the visco-resistive MHD 
simulation. The format is similar to that used in Tables [UlSl Column (1) specihes the model 
applied. The hrst letter V denotes “visco-resistive.” The second to fourth numbers denote the 
logarithmic of the plasma beta, the viscosity (multiplied by 10^), and the magnetic Prandtl 
number, respectively. Columns (2) and (3) list the viscosity and resistivity, respectively. 

Figure [10] shows the correlation of statistical averages in the turbulent state. The 
Maxwell stress is smaller than the magnetic energy by a factor of 2 — 3, and larger than 
the Reynolds stress by a factor of 5 — 10. The correlation is consistent with the ideal MHD 
case (Figure [2]). The magnetic held is anisotropic, and it is preferentially enhanced along 
the ^/-direction. 

Figure [TT] presents the statistical average of the Maxwell stress over six runs: V7-0-1, V?- 
1-1, V7-2-1, V7-2-2, V7-2-4, and V7-1-0. In Figure fTTf ah the Maxwell stress is very weakly 
dependent on the pressure when the magnetic Prandtl number is hxed. Note that the ideal 
MHD simulation with the same numerical scheme shows positive correlation (cross symbols 
in Figure [6|(b) and red-colored symbols in Figure [9|), wherein the numerical magnetic Prandtl 
number is an increasing function of the pressure. We conhrm that the visco-resistive MHD 
simulation with a different Riemann solver yields a consistent result. A similar tendency 
is also observed in the case of the ideal MHD simulation with the F1-HLL-F2-HLL scheme 
(Fignre[6](a) and blue-colored symbols in Figure [9]), wherein the numerical magnetic Prandtl 
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number is nearly constant. Therefore, we conclude that the Maxwell stress in the MRI- 
induced turbulence at a constant magnetic Prandtl number is almost pressure-independent 
within the explored parameter range. Figure [rTT bl shows the Maxwell stress as a function of 
the magnetic Prandtl number. The stress curve roughly obeys the expression wm oc 
T he power is close to the correlation obtained in Section ITTI This consistency indicates that 
the estimation of our numerical magnetic Prandtl number and hence the conclusion in Section 
13.II are quantitatively valid. 

In Table El we hnd an unexpectedly high level for the Maxwell stress especially at the 
cases with high viscosity (e.g., V6-4-4). Its origin is discussed in Appendix [Bl 


4. DISCUSSION 


4.1. Magnetic Prandtl Number Dependence 


We discuss the magnetic Prandtl numbe r dependence of the s atura tion level of the MRI 
from a physical point of view. As is shown by iLesnr fc Longarettil (120071 ) and IPessah fc Chan 
(j2008|), the linear growth rate of the axisymmetric incompressible MRI is an increasing func¬ 
tion of the Reynolds and magnetic Reynolds numbers, thus not of the magnetic Prandtl num¬ 
ber. Consequently, we consider dissipation by secondary instabilities as a possible candidate. 
It has been widely recognized that the so-called channel flow of the primar y MRI mode be¬ 
come s unstable against secondary three-dimensional parasitic instabilities (ICoodman fc Xu 
1994ji. The instabi l ities are related to the Kelvin-Helmholtz (KH) and tearing modes. 


Pessah fc CoodmanI (120091) have carried out a linear analysis of the secondary parasitic in¬ 


stability under the conhguration of the primary MRI mode and have shown the relative 
contribution of KH and tearing modes to the parasitic instability. Although they present a 
rigorous description for the secondary modes, we simply consider the KH and tearing insta¬ 
bilities independently to derive their Reynolds (Rg) and magnetic Reynolds (Rm) number 
dependence. Details of the linear analysis of the KH and tearing instabilities are described 
in Appendix O 


Figure IT^ al shows the linear growth rate of the KH instability. We assume that the 
Alfven velocity va is 1.5 times larger than the half velocity jump and an angle between 
the shear flow and magnetic held 9 = 60°. The linear growth rate of the KH instability is 
almost independent of the Reynolds and magnetic Reynolds numbers when Re,Rm > 10^. 
A similar conclusion can be drawn for different combinations of va and 6. Therefore, the 
KH instability does not appear to cause the magnetic Prandtl number dependence. 


On the other hand, the tearing instability shows magnetic Prandtl number dependence 
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in Figure IT^ b.c). The linear growth rate is a decreasing function of the magnetic Reynolds 
nu mber iFigure flW bR and an increasing function of the Reynolds number (Figure fl^ cl) 
(cf., |Porcellilll987|) . The Alfven velocity with respect to the dominant magnetic held By is of 
the order of 10“^, thereby yielding ~ = 10^ during the saturation state of the MRI. 

The growth rate can be htted by ytear oc Rm^'^R^fJ^ = Rrm^^Rm^"^ around Re ~ Rm = 10^, 
a decreasing function of the magnetic Prandtl number. 


The saturation state is asse ssed by equating the growth rate of the primary (MRI) and 
the secondary (tearing) modes (jPessah fc Goodman! 1200911 . In the course of the MRI, the 
dominant magnetic held By stretched along the y-direction will become unstable against the 
tearing mode. The condition for the saturation may be 


Ttear OC 


MRI 




m m 


7mri- 


Using the relation 7mri/^mri = vaz = constant (e.g., iMasada fc Sanoll2008l) . we obtain 


B,'‘ a; 


(16) 


(17) 


Since the Maxwell stress is proportional to the magnetic energy (Figure flU]) . equation (ITT)) 
supports the simulation result of wm oc Pr^ 


.5-1.0 

m 


The linear analysis of equation flTTl) indicates that the saturation level depends on the 
magnetic Reynolds number as well as the magnetic Prandtl number. However, the de¬ 
pendence on the magnetic Reynolds number is less pronounced than that on the magnetic 
Prandtl number within the explored parameter range in our visco-resistive simulation. The 
mag netic Reynolds numb er d ependence of t he sa turation level has been investigated by 




Fleming et al.l ( 2000l) and lSano fc Stonel (120021) by means of resistive MHD simulations 


Sano fc Stonel (120021) have shown that the saturation level of the Maxwell stress decreases 
when v\/riVL < 1, and further, the level is nearly constant when v\/7]VL > 1. The magnetic 
Reynolds number dependence does not appear to obey a single power-law relation. 

The rigorous linear analysis of the secondary parasitic instability under the primary MRI 
mode has shown that the fastest parasitic instability is associat ed with the KH mode, no t 
the tearing mode, especially at high magnetic Prandtl number (jPessah fc Goodma.nll2009l) . 
By equating the growth rates of the primary MRI and the secondary parasitic instability, 
the saturation level is almost independent of the magnetic Prandtl number when the spatial 
domain is unlimited to permit the fastest primary and secondary instabilities. However, they 
have shown that the saturation level weakly depends on the magnetic Prandtl number under 
the domain size of {L^.^ L^, L^) = 1 x 4 x 1. The growth of the KH instability is suppressed 
by limiting the domain in the horizontal direction whereas the domain is sufficient in the 
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azimuthal direction so as to permit the tearing instability. This situation may overestimate 
the relative contribution of the tearing mode to the parasitic instability, thereby causing the 
magnetic Prandtl number dependence. In this regard, our interpretati on of the simulation 
based on the simplihed linear analysis follows the rigorous analysis by iPessah fc Goodman 
(I 2 OO 9 I) . As they have pointed out, it is necessary to use a sufficiently wide simulation domain 


(at least, Ly > 2Lz) to accommodate the fastest parasitic modes. The (numerical) 
magnetic Prandtl number dependence of the saturation level may be moderated provided 
we use a w ide domain i n the horizontal direction so as to permit the fastest parasitic (KH) 
instability. iBodo et al.l fj2008|) have investigated the aspect ratio dependence of the saturation 
level of the MRl in the unstratihed shearing box. The channel flow that is well developed in 
the case with an aspect ratio of unity L^/L^ = 1 disappears at an aspect ratio larger than 
unity, with which the KH instability is expected to disrupt the flow. Insufficient domain 
size in the horizontal direction may lead to overestimating the saturation level of the MRl- 
induced turbulence. 

The correlation between the magnetic Prandtl number and Maxwell stress as obtained 
bv lhesur fc Longarettil (120071) {wm oc is weaker than our result {wm oc 

indicating that other parameters affect th e dependence. This d iscrep ancy may be due to the 
different range of the Reynolds number. Ihesur fc Longarettil (120071) explored the magnetic 
Prandtl number dependence in a wide range of the Reynolds number = 200 — 6400 by 
virtue of a pseudo-spectral incompressible code whereas our si r nulati ons are Re < 2000. The 
magnetic Prandtl number dependence in Ihesur Sz Longarettil (120071) is more prominent for 
lower Reyn olds number (Fig. 1 0 in th eir paper). A similar tendency is found in the linear 
analysis bv IPessah fc GoodmaiJ (120091) . although the dependence is less pronounced than in 
the case of nonlinear simulation results. 


4.2. Gas Pressure Dependence 

The pressure-independent saturation level of the MRl is indicative of the incompressible 
turbulence in the intermediate to high pressure range. On the other hand, the saturation 
level shows positive correlation with the gas pressure at very low pressures (the leftmost 
symbols in Figures l5](a), IHl and 19]) regardless of the choice of the scheme. This may 
indicate the role of physical dissipation for the saturation of the MRl, which is irrelevant 
to numerical dissipation scales. The magnetic pressure is amplihed by a factor of ~ 100 in 
our ideal MHD simulations regardless of the initial gas pressure. Therefore, the plasma beta 
at the saturation is close to or smaller than unity in the lowest pressure case (/3 = 10^). 
The flow speed of ~ va driven by the MRl becomes supersonic with f3 < 1. Consequently, 
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compressibility may play an important role in the saturation of the MRI. For example, slow¬ 
mode shocks contribute to the dissipation of the current and eventually lead to the decrease 
of the Maxwell stress with decreasing the pressure (increasing compressibility). As can be 
found in Figure the time variation of the stress is signihcantly suppressed for the lowest 
pressure case. This implies that a large-scale channel flow is inhibited from growing in large 
amplitudes by compressibility. Subsequently, it will inhibit the secondary parasitic instability 
because its growth rate is proportional to the amplitude of the primary mode. Alternatively, 
Latter et al.l (1200911 have argued that non-linear multi-mode interactions by turbulent mixing, 
not the linear parasitic instability, is responsible for the saturation of the MRI in the case 
when the primary mode cannot reach large enough amplitudes. Although the compressible 
case (/5 = 10^) in our simulation is not directly applicable to actual accretion disks because 
the vertical size is considerably larger than the pressure scale height, L ^ \/‘^P/p/^t the 
simulation of the MRI with relatively low plasma beta is of interest to the fundamental 
study of compressible MHD turbulence. The saturation mechanism of the compressible 
MRI-induced turbulence will be investigated in detail in future. 

Although the saturation level of the MRI is found to be almost pressure-independent 
when the magnetic Prandtl number is constant, it does not necessarily hold true in actual 
plasma environments. The saturation level may depend on the pressure because the magnetic 
Prandtl number can change with density, temperature, and hence pressure. For example, 
the magnetic Prandtl number in weakly magnetized, collisional plasmas is estimated to be 
Pvrr, ~ 10 ~^T ^ /n for fully ioniz e d cas es and Pr^ ~ lO^T^/n for partially ionized cases 


flSpitzeiill962l) . iPotter fc BalbusI (1201411 have considered the magnetic Prandtl number for 
radiation-dominated disks, Pr^ 10iOTiiA/(fi;n2), where k = (0.4-|-8nT cm^g ^ is 
the opacity for bound-free absorption in a partially ionized gas. 


SUMMARY 


Magnetohydrodynamic (MHD) turbulence has been extensively studied via ideal MHD 
simulations. In general, ideal MHD simulations of turbulence with various schemes are not 
identical because viscous and resistive dissipation scales differ among schemes. We have 
investigated the saturation level of the MHD turbulence induced by the magnetorotational 
instability (MRI) by means of the ideal MHD simulation. We use Godunov-type schemes 
with various reconstruction functions and Riemann solvers in order to assess effects of nu¬ 
merical dissipation. Upon revis iting the ga s press ure dependence of the saturation level of 
the Maxwell stress as studies bv ISano et ahl (120041) . we have shown that the saturation level 
positively correlates with the pressure only at very low pressure, in which the MRI-driven 















19 


flow becomes compressible. Except for this case, we have failed to evaluate the gas pres¬ 
sure dependence of the saturation level because numerical viscous and resistive scales vary 
according to the pressure as well as the numerical scheme. 


We have estimated the numerical magnetic Prandtl number Pr^ (=kinematic viscosity 
///resistivity 77 ), and we have shown that the saturation level is more sensitive to the numerical 
magnetic Prandtl number than the pressure. Since the numerical magnetic Prandtl number 
itself depends on the pressure, the pressure dependence of the saturation level is indistin¬ 
guishable from the numerical magnetic Prandtl number dep e ndenc e. Therefore, we conclude 
that the gas pressure dependence obtained by ISano et al.l (120^ could be a consequence 
of the pressure-dependent numerical magnetic Prandtl numbeio. We have also conducted 
the same runs with the ZEUS code, and we have conhrmed that both the saturation level 
and the numerical magnetic Prandtl number are almost pressure-independent (Appendix 
El). One should recognize the numerical magnetic Prandtl number of the applied numerical 
code and its impact on the result. Also, the result strongly suggests the need for the ex¬ 
plicit use of physical viscosity and resistivity to control the m agnetic Prandtl number f or the 
MR I-induced turb ul ence simulation as poin t ed ou t by e.g., iFromang fc Papaloizoul (120071) 
and ISilversI (120081 ). iFromang fc Papaloizoul (120071) showed that solutions of the ideal MHD 
simulation of the MRI-in duced turbulence with zero net magnetic flux depend on numerical 
resolution. ISilversI (120081) showed the numerical resolution dependence of the saturation level 
even with hnite net flux cases, in which the minimum scale is limited by the critical wave¬ 
length of the MRI with respect to the hnite ambient magnetic held, A = 27rn^/y^f2 oc Bq. 
In addition to them, we have revealed that the saturation level depends on the numerical 
magnetic Prandtl number, which is more sensitive to characteristics of the applied scheme 
rather than numerical resolution. 


Consequently, we have carried out a visco-resistive MHD simulation to investigate the 
parameter dependence of the saturation level of the MRI-induced turbulence. The satu¬ 
ration level of the Maxwell stress depends on the physical magnetic Prandtl number. It 
is almost independent of the gas pressure when the magnetic Prandtl number is constant 
and compressibility is weak. The positive correlation between the Maxwell stress and the 
magnetic Prandtl number may be an indication of the saturation mechanism. We consider 
from linear analysis that the secondary tearing instability saturates the growth of the MRI, 
and consequently, it causes the magnetic Prandtl number dependence. 


^In fact, the numerical magnetic Prandtl number for the scheme used in their paper is found to be an 
increasing function of the pressure. 
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A. IMPLEMENTATION OF THE MULTI-STATE RIEMANN SOLVER TO 
THE UPWIND CONSTRAINED TRANSPORT METHOD 


We consider the two-dimensional ideal MHD equations, 


dU dF dG 

dt ^ dx ^ dy 

(Al) 

dBx dEz 
dt dy 

(A2) 

dBy dEz _ 
dt dx 

(A3) 


where U = (p, Ux, Uy, u^, e, Bz) denotes conservative variables (e the total energy), F(U, Bx, By) 
and G{U, Bx, By) the corresponding fluxes along x and y directions, and Ez = UyBx —UxBy 
the electric held, respectively. The extension to three dimension is straightforward. The 
equations are discretized into a hnite volume formulation as 


dUij _ ^1+1/2^ ~ Fi_i/2J Gjj+i/2 — GiJ-l/2 

dt Ax Ay 

dBx-i-l/2,j _ Ez-i-l/2,j+l/2 — Ez-i-ll2,j-ll2 

dt Ay ’ 

dBy.ij_l/2 ^ Ez-i+ll2,j-ll2 - Ez-i-ll2,j-ll2 
dt Ax ’ 


(A4) 

(AS) 

(A6) 


where Uij dehned at cell centers are averaged over the area and Bx-i-ii 2 ,j, By.ij_i /2 dehned at 
cel l faces are averaged alon g the orthogonal line. This is the staggered grid system employed 
by Evans fc Hawlev fjl988 h 


We solve the above equations with an upwind scheme based on Riemann solvers. First, 
we interpolate the in-plane magnetic helds to cell centers as cell-averaged representations. 
For example, the second, fourth, and sixth order approximations are respectively expressed 
as 

Ey,i,j = 2 A Ey-iJ+l/2) , 

^ {-By-ij-3l2 + 13Ry;iJ-l/2 A 1 J+ 1/2 - By.ij+ 3 / 2 ) , 


(A7) 

(AS) 
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B. 




(lli?y;ij_5/2 - 93i?y;ij_3/2 + 802i?y.jj_i/2 
+802i?j^.jj + l/2 — 93i?j^;jj+3/2 + ll-Bj/;ij+5/2) , 


(A9) 


where we assume uniform grid spacing. We adopt the sixth order approximation. Since all 
the variables Wi^j = (Uij, By-ij) are prepared at cell centers, they are interpolated 
along the x-direction by an arbitrary reconstruction function Ti with the degree of 2s + 1 
to obtain the left and right states at right and left faces, 




Lx 


j, 


*+1/2 J’ 


i+s,j) 


(AlO) 


where the superscripts Lx, Rx indicate the left and right states along the x-direction. Note 
that the interpolation of Bj. (normal component) to cell faces is not required because it is 
already dehned there. Subsequently, we solve the Riemann problem at cell faces with an 
arbitrary one-dimensional Riemann solver TZi to determine the variable W 1 ^ 1 / 2 ,j as well as 
the upwind numerical flux R’j_i/ 2 j- For example , the HLLD approximate Riemann solver 
determines them to satisfy F’hlld = -FChFHLLD) (IMivoshi fc Kusanol 120051) . 


The above procedure is the one-dimensional reconstruction along the x-direction at 
Hj. To determine the electric field Ez-i-i/ 2 ,j-i /2 at cell edges (numerical flux for the in-plane 
magnetic fields), we perform the same procedure along the y-direction at Xj_i/ 2 - The variable 
Wi - 1 / 2 ,j is interpolated along the |/-direction by a reconstruction function R 2 to obtain their 
left and right states at right and left edges. 


i-l/2,j+l/2 


’ ** i- 


Ry 

i-l/2j-l/2 


•(— J ^2 (FFj_l/ 2 j-S 5 • • • 5 FFi-l/ 2 j 5 ■ ■ ■ 5 FFj_i/ 2 j+s) • (All) 


Subsequently, we solve the Riemann problem at cell edges with a Riemann solver IZ 2 to 
determine the upwind electric field. 


Ez-i-ll2,j-ll2 ^ R 2 



Ly 

i-l/2j-l/2 


,w 


Ry 

i-l/2j-l/2 


Here, a multi-state Riemann solver is available for 7^2- 


(A12) 


For symmetry, we perform the same procedure that involves the reconstruction along 
the |/-direction at x* to obtain VFjj_i /2 followed by the reconstruction along the x-direction 
at l/j- 1/2 to determine 


W j+i/ 2 j-l/ 2 ) i_i/ 2 J-l /2 ^ L 2 (vr i-sj- 1 / 2 , • • • , 

Ez-i-l/2,j-l/2 ^ R 2 (FF^i"-l/2,i-l/2) FFf!)/2,j-l/2) • 


(AM) 


Finally, we use the arithmetic average of Equations flA12D and flAMp as the numerical flux 
in Equations flASp and flA6D . A similar method has been recently prop osed to advance the 
Maxwell equation in kinetic plasma simulations (IMinoshima et al.l 120151 ). 
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Our method is not identical t o the upwind constrained transport method originally pro¬ 
posed by iLondrillo fc del Zannal fl2004l) in that they use a two-dimensional HLL Riemann 
solver whereas we use an arbitrary one-dimensional Riemann solver dimension-by-dimension. 
We consider that the essence is common between the two methods; calculation of the electric 
held at cell edges to retain the multi-dimensional upwind property. Our method can easily 
incorporate multi-state Riemann solvers by virtue of successive one-dimensional reconstruc¬ 
tions. 


The accuracy of the magnetic held is determined by J^ 2 i'T ^2 whereas that of the huid is 
by In this paper, we use various combinations of TZi and ”^2 to assess ehects 

of numerical dissipation in the ideal MHD simulation of the MRI-induced turbulence. 


B. TRANSITION BETWEEN LAMINAR AND TURBULENT STATES 


As can be found in Table El the Maxwell stress is unexpectedly high at the cases with 
high viscosity. For example, the Maxwell stress is nearly constant for the range of u = rj = 
(0.5 — 2) X 10“^ (V?-0-l,V?-l-l,V?-2-l), but it is suddenly increased at u = rj = 4 x 10“^ 
(V7-4-1) and is fairly sensitive to the gas pressure. Figure 113] shows the time prohle of the 
Maxwell stress and the vertical magnetic held energy R^/Stt normalized by the pressure in 
the run denote d as V4-4-1. There are m any quasi-periodic bursts with exponential growth 
as reported bv iLesur fc Longarettil (1200711 . Obviously, time averaging is meaningless in such 
a situation. During the burst growth, the vertical magnetic held energy is nearly equal to 
the initial level and the growth rates of bursts are almost equal. Therefore, the bursts are 
expected to be essentially identical to the initially unstable mode. The periodic appearance 
of the burst suggests that the solution is laminar rather than turbulent. 

The origin of the bu rst is understood from linear theory. The dispersion relation (e.g., 
Lesur fc Longarettil 1200711 is analytically solved at = 1 as 


7 

n 


1 f ^ 

V D 


+ 



+ {2-qf- 


kVA 

n 


- (2 - g) 


1/2 


(Bl) 


where S = v\/rjVL{= v^/uQ) denotes the Lundquist number. The growth rate and the 
wavenumber of the fastest growing mode (FGM) are approximated as 


vD/fgm 2 \ 2S J \ D / FGM l + 2g/S' 4 


(B2) 
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for S'1, and 


12/PGM 


2 


2q 


at 


kVA\ _ S 

^ / PGM 


(B3) 


for S' 1. The critical wavenumber for the instability is kcritVA/^ = \/2g for S' S> 1 (ideal 
limit) and kaitVA/^ = — <?)S' for S' 1 . 

The instability condition is characterized by a combination of the Lundquist number S 
and the minimum wavenumber limited by the vertical size k^i^ = 27r/L. When the minimum 
wavenumber is considerably smaller than that of the FGM, the system involves multiple 
unstable modes, and subsequently, the turbulence is successfully sustained by the MRI. On 
the other hand, only the longest wavelength mode grows when the minimum wavenumber is 
larger than that of the FGM and smaller than the critical wavenumber. In such a situation, 
we o bserve only the largest laminar flow that is marginally unstable flLesur fc Longaretti 
2007 1. In the runs at = 1 (V?-0-l,V?-l-l,V?-2-l,V?-4-l), the condition for the laminar 
regime for F <C 1 , 


1 kmirJ^ 


V2 


< 


va 


< 


<1 


Q 


= 73, 


(B4) 


is satished only with z/ = 4 x 10“^ (V7-4-1), in which the quasi-periodic burst with k = fcmin 
is actually observed. The growth rate of 7 = O.IG is in good agreement with the exponential 
growth of bursts. 

The laminar regime also appears in the high magnetic Prandtl number cases (e.g., 
runs V7-4-2 and V7-4-4). This is expected as per linear theory in the viscous limit (z/ S> 
v\/VL,r] = 0). The wavenumber of the fastest growing mode decreases with the viscosity, 
fcpGM = \/ 2{2 — q)Vt/ y ^ but the critical wave number remains unchanged, fccrit = \/ 2 gG/nA 
fjPessah fc GhanI l2008l: iMasada fc Sand l2008l ). Therefore, the largest flow is unstable re¬ 


gardless of the viscosity. The condition for the laminar regime in the viscous limit is given 
by 


Va 


L 

< — < 


271 Y 2(2 - q)n 


V 


(B5) 


The exponential burst is most likely to be observed at high viscosity. This is not the 
case in the resistive limit (z/ = 0,7 ^ v\/VL) because the critical waven umber as well as 
the w avenumber of the fastest growing mode decrease with the resistivity fjSano fc Mivama 
1999h . 
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C. LINEAR ANALYSIS OF VISCO-RESISTIVE KELVIN-HELMHOLTZ 

AND TEARING INSTABILITIES 


In the incompress ible visco-resistive MHD, th e linearized equati ons of the KH instability 
f) Chandrasekhar! 1 196 ll) and the tearing instability flFurth et al.lll963l) are respectively written 


as 


7 khA = ik [cos 6iIj — FA] + R^AA, 
7khA'0 = ik F"'i(j — FA'l/j + j cos 


BAA 


+ R-^AA^/j, 


(Cl) 


and 


r 7tear^ = Fi/j + RjAA, 

{ 7tearAV^ = k^ {F"A - FAA) + R-^AAij, 


(C2) 


where A and denote the flux and stream functions, Re and Rm the Reynolds and magnetic 
Reynolds numbers, A = —k^ + (F/dz^, F = tanh(z) the prohle of the background shear 
flow (KH) or magnetic held (tearing), and F” = d?F/dz^. The parameters va and vq 
denote the Alfven velocity and the half velocity jump, 9 denotes the angle between the 
shear how and magnetic held, and we assume that the KH wave is parallel to the shear 
how. The in-plane magnetic held is included in the KH instability, and the background 
how is ignored in the tearing instability for simplicity. In equation flCip . the variables are 
normalized by uq and the width of the shear layer <5. The Reynolds and magnetic Reynolds 
numbers are expressed as Re = VoS/u and Rm = Vod/rj, respectively. In equation flC2D . the 
variables are normalized by the Alfven velocity and the width of the current layer 6, thus the 
Reynolds and magnetic Reynolds numbers are expressed as Re = VA^jv and Rm = vaS/t], 
respectively. We numerically obtain the eigenvalues and eigenfunctions with the boundary 
condition AA = Aip = 0 at \z\ = ±10. 


D. SUPPLEMENTARY SIMULATIONS USING THE ZEUS CODE 


Using the ZEUS code ( Stone fc Norman 1992alji3 ). we ran a set of ideal MHD simulations 
to investigate the gas pressure dependence of the saturation level of MRI. The parameters 
are the same as those in the text except for the net vertical held Bq, which is larger by a 
factor of \/4^. (There is no particular reason for this choice of Bq.) Figure [131(a) shows 
the Maxwell stress (averaged from 20 to 100 orbits) vs. gas pressure. The Maxwell stress is 
almost independent of gas pressure except the smallest pressure case. Then we evaluated the 
numerical magnetic Prandtl numbers for these simulations in the same manner as explained 
in the text except that each number is calculated from a statistical average of four simulations 
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restarting from different epochs (20, 40, 60, and 80 orbits). Figure fldlf bl shows that the 
numerical magnetic Prandtl number is roughly constant regardless of gas pressure. These 
results support the statement in the text that the saturation level does not depend on gas 
pressure (except for very low plasma beta cases) when (numerical) magnetic Prandtl number 
is constant. 
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Grid size 


Fig. 1.— Li error norm of the two-dimensional Alfven wave propagation problem as a 
function of the grid size. The diamonds and triangles are obtained with the third and fourth 
order Runge-Kutta time integration, respectively. The dotted, dashed, and dash-dotted lines 
represent the third, fourth, and hfth order of accuracy, respectively. 
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Magnetic energy Reynolds stress 


Fig. 2.— Correlation of statistical averages in ideal MHD simnlations of the MRI-induced 
tnrbulence between (a) the magnetic energy and the Maxwell stress, and (b) the Reynolds 
stress and the Maxwell stress. 
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Fig. 3.— Time profile of the Maxwell (black line) and Reynolds (gray line) stresses in the 
runs (a) I2-WRZH, (b) I4-WRZH, and (c) I6-WRZH. 
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Fig. 4.— Spatial profile of the |/-component of the magnetic held at (a) peak and (b) decay 
periods in the run I4-WRZH. Left panels show the x — z plane ai y = —2, and right panels 
show the y — z plane at x = 0.5. White arrows indicate the direction of how (background 
shear is subtracted). 
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Fig. 5.— (a) Statistical average of the Maxwell stress as a function of the gas pressure. The 
HLLR and HLL Riemann solvers are used for the fluid and magnetic held updates, respec¬ 
tively. The various symbols represent simulation results obtained with different reconstruc- 
tion functions. T he dash-dotted line represents the correlation wm oc as obtained by 


Sano et ahl (120041) . (b) Time prohle of the Maxwell stress in the runs I6-WRWH (black line) 


and I6-WRZH (gray line), corresponding the cross and asterisk, respectively, at P = 312.5 
in panel (a). 
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(a)HLL+HLL 


(b)HLLR+HLLR 
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Fig. 6.— Statistical average of the Maxwell stress as a function of the gas pressure. The 
(a) HLL and (b) HLLR Riemann solver are used. The various symbols represent simulation 
results obtained with different reconstruction functions. The dash-dotted line represents the 
correlation wm oc as obtained by 


Sano et al.l te004j) . 
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Fig. 7.— Example of time profile of the numerical magnetic Prandtl number. The WENO- 
HLLR-WENO-HLL scheme is used. The various lines represent simulation results obtained 
with different initial plasma beta values. 
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Fig. 8.— Statistical average of the Maxwell stress (a) as a function of the numerical magnetic 
Prandtl number, and (b) as a function of the gas pressure and the numerical magnetic Prandtl 
number 
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Fig. 9.— Statistical average of the Maxwell stress as a function of the gas pressure. Open and 
hlled diamonds correspond to low and high resolution runs, respectively. The various colors 
represent simulation results obtained with different combinations of the R iemann solv e r. Th e 
dash-dotted line represents the correlation wm oc as obtained bv ISano et all (120041) . 
Note that the symbol for the low resolution with the WENO-HLLR-WENO-HLLR scheme 
(red) is overlapped with that for the high resolution at P = 31.25. 
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Magnetic energy Reynolds stress 




Fig. 10.— Correlation of statistical averages in visco-resistive MHD simulations of the MRI- 
induced turbulence between (a) the magnetic energy and the Maxwell stress, (b) the Reynolds 
stress and the Maxwell stress, (c) B‘1 and and (d) B‘1 and By. 
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Fig. 11.— Statistical average of the Maxwell stress in visco-resistive MHD simulations as a 
function of (a) the gas pressure and (b) the physical magnetic Prandtl number. The various 
symbols represent simulation results obtained with different combinations of the viscosity 
and resistivity. 
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Fig. 12.— Linear growth rate in the incompressible visco-resistive MHD. (a) The Kelvin- 
Helmholtz instability as a function of the Reynolds number. (b,c) The tearing instability as 
a function of the magnetic Reynolds number and the Reynolds number. 
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Fig. 13.— Time profile of the Maxwell stress (red line) and the vertical magnetic held energy 
(blue line) in the run V4-4-1. 
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Fig. 14.— Pressure dependence of the saturation level of MRI (a) and the numerical magnetic 
Prandtl number (b) based on ideal MHD simulations using the ZEUS code. 
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Table 1: Ideal MHD simulations of the MRI. The HLLR and HLL Riemann solvers are used 


for the fluid and magnetic held updates, respectively. 


(1) Model 

(2) FI 

(3) F2 

{4)/3 

(5) a, 

(6) Em X 10“ 

(7) Ek X 10“ 

(8) wm X 10^ 

(9) We X 10“ 

(10) a X 10“ 

(11) Pr^,num 

F1-HLLR-F2-HLL 

I2-WRWH 

WENO 

WENO 

10“ 

3.125 X 10-“ 

4.05 

1.23 

1.46 

0.357 

58.1 

0.76 

I3-WRWH 

WENO 

WENO 

10^ 

3.125 X 10“^ 

6.49 

2.60 

2.76 

0.745 

11.2 

0.43 

I4-WRWH 

WENO 

WENO 

10-* 

3.125 X 10" 

5.23 

2.43 

2.32 

0.642 

0.948 

0.51 

I5-WRWH 

WENO 

WENO 

10“ 

3.125 X 10^ 

4.19 

1.72 

1.81 

0.411 

0.0710 

0.45 

I6-WRWH 

WENO 

WENO 

10“ 

3.125 X 10“ 

2.93 

1.03 

1.07 

0.258 

0.00426 

0.36 

I2-WRZH 

WENO 

WZMP 

10^ 

3.125 X 10-“ 

4.82 

1.42 

1.75 

0.438 

70.3 

1.6 

I3-WRZH 

WENO 

WZMP 

10^ 

3.125 X 10“^ 

8.07 

3.21 

3.52 

0.946 

14.4 

0.88 

I4-WRZH 

WENO 

WZMP 

10“ 

3.125 X 10" 

7.63 

3.25 

3.47 

0.871 

1.39 

1.4 

I5-WRZH 

WENO 

WZMP 

10“ 

3.125 X 10^ 

9.28 

3.61 

4.31 

0.957 

0.169 

1.4 

I6-WRZH 

WENO 

WZMP 

10“ 

3.125 X 10“ 

12.5 

4.30 

5.53 

1.21 

0.0215 

0.84 

I3-ZRZH 

WZMP 

WZMP 

10^ 

3.125 X 10“^ 

7.13 

3.05 

3.07 

0.851 

12.7 

0.51 

I4-ZRZH 

WZMP 

WZMP 

10“ 

3.125 X 10" 

5.28 

2.77 

2.45 

0.684 

1.00 

0.82 

I5-ZRZH 

WZMP 

WZMP 

10“ 

3.125 X 10^ 

5.38 

2.52 

2.47 

0.580 

0.0977 

0.74 

I6-ZRZH 

WZMP 

WZMP 

10“ 

3.125 X 10“ 

4.32 

2.04 

1.95 

0.493 

0.00782 

0.48 
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Table 2: Ideal MHD simulations of the MRI. The HLL Riemann solver is used for both the 


fluid and magnetic held updates. 


(1) Model 

(2) FI 

(3) F2 

(4)/J 

(5)fo 

(6) Em X 10" 

(7) Ek X 10" 

(8) wm X 10^ 

(9) WR X 10" 

(10) Q X 10" 

(11) Pr^,num 

Fl-HLi:^F2-HLL 

I2-WHWH 

WENO 

WENO 

10" 

3.125 X 10-" 

2.01 

1.04 

0.848 

0.307 

37.0 

0.95 

I3-WHWH 

WENO 

WENO 

10^ 

3.125 X 10“^ 

6.55 

2.26 

2.79 

0.655 

11.0 

0.91 

I4-WHWH 

WENO 

WENO 

10^ 

3.125 X 10“ 

7.58 

2.84 

3.50 

0.781 

1.37 

0.86 

I5-WHWH 

WENO 

WENO 

10“ 

3.125 X 10^ 

6.93 

1.96 

2.97 

0.496 

0.111 

0.94 

I6-WHWH 

WENO 

WENO 

10“ 

3.125 X 10" 

5.78 

1.31 

2.20 

0.359 

0.00819 

0.92 

I2-WHZH 

WENO 

WZMP 

10^ 

3.125 X 10“" 

3.19 

1.15 

1.27 

0.349 

51.8 

1.3 

I3-WHZH 

WENO 

WZMP 

10^ 

3.125 X 10“^ 

8.24 

2.77 

3.58 

0.846 

14.2 

1.2 

I4-WHZH 

WENO 

WZMP 

10^ 

3.125 X 10“ 

9.02 

3.04 

4.16 

0.877 

1.61 

1.4 

I5-WHZH 

WENO 

WZMP 

10“ 

3.125 X 10^ 

10.8 

3.03 

5.06 

0.832 

0.189 

1.3 

I6-WHZH 

WENO 

WZMP 

10“ 

3.125 X 10" 

11.7 

2.59 

5.07 

0.782 

0.0187 

1.6 

I3-ZHZH 

WZMP 

WZMP 

10^ 

3.125 X 10“^ 

8.43 

3.22 

3.66 

0.922 

14.7 

0.90 

I4-ZHZH 

WZMP 

WZMP 

10^ 

3.125 X 10“ 

10.0 

3.98 

4.61 

1.01 

1.80 

1.1 

I5-ZHZH 

WZMP 

WZMP 

10“ 

3.125 X 10^ 

6.30 

2.54 

3.02 

0.579 

0.115 

1.3 

I6-ZHZH 

WZMP 

WZMP 

10“ 

3.125 X 10" 

6.60 

2.43 

3.10 

0.582 

0.0118 

0.99 
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Table 3: Ideal MHD simulations of the MRI. The HLLR Riemann solver is used for both the 


fluid and magnetic held updates. 


(1) Model 

(2) FI 

(3) F2 

(4)/J 

(5) a 

(6) Em X 10" 

(7) Ek X 10" 

(8) wm X 10^ 

(9) wj, X 10" 

(10) a X 10" 

(11) Pr^,num 

F1-HLLR-F2-HLLR 

I2-WRWR 

WENO 

WENO 

10" 

3.125 X 10-" 

6.32 

1.68 

2.27 

0.535 

90.4 

2.7 

I3-WRWR 

WENO 

WENO 

10^ 

3.125 X 10“^ 

11.4 

3.87 

4.93 

1.19 

19.7 

2.6 

I4-WRWR 

WENO 

WENO 

10" 

3.125 X 10" 

14.3 

4.90 

6.54 

1.40 

2.54 

3.2 

I5-WRWR 

WENO 

WENO 

10® 

3.125 X 10^ 

27.6 

7.38 

12.0 

1.97 

0.447 

4.3 

I6-WRWR 

WENO 

WENO 

10“ 

3.125 X 10" 

196 

31.0 

70.4 

8.95 

0.253 

4.9 

I2-WRZR 

WENO 

WZMP 

10^ 

3.125 X 10“" 

6.60 

1.73 

2.34 

0.553 

93.4 

3.6 

I3-WRZR 

WENO 

WZMP 

10^ 

3.125 X 10“^ 

14.3 

4.60 

6.11 

1.41 

24.2 

2.9 

I4-WRZR 

WENO 

WZMP 

10" 

3.125 X 10" 

21.3 

6.31 

9.11 

1.76 

3.47 

4.1 

I5-WRZR 

WENO 

WZMP 

10“ 

3.125 X 10^ 

42.7 

10.5 

17.8 

2.79 

0.659 

6.7 

I6-WRZR 

WENO 

WZMP 

10“ 

3.125 X 10" 

354 

53.9 

126 

17.2 

0.457 

6.8 

I3-ZRZR 

WZMP 

WZMP 

10^ 

3.125 X 10“^ 

16.0 

5.49 

6.82 

1.59 

27.3 

2.8 

I4-ZRZR 

WZMP 

WZMP 

10" 

3.125 X 10" 

22.3 

7.13 

9.66 

1.81 

3.69 

3.3 

I5-ZRZR 

WZMP 

WZMP 

10“ 

3.125 X 10^ 

21.7 

6.58 

9.71 

1.63 

0.364 

5.6 

I6-ZRZR 

WZMP 

WZMP 

10“ 

3.125 X 10" 

53.6 

12.6 

21.9 

3.26 

0.0806 

6.9 






-45 - 


Table 4: High resolution ideal MHD simulations of the MRI. 


(!) Model 

(2) R1 

(3) R2 

WP 

(5) P„ 

(6) X 10" 

(7) Ek X 10" 

(8) WM X 10^ 

(9) WR X 10^ 

(10) a X 10" 

(11) Prm,num 

H2-WRWH 

HLLR 

HLL 

10" 

3.125 X 10-" 

5.27 

1.65 

1.87 

0.463 

74.9 

0.53 

H3-WEWH 

HLLR 

HLL 

10^ 

3.125 X 10“^ 

8.77 

3.46 

3.69 

0.913 

14.8 

0.41 

H4-WRWH 

HLLR 

HLL 

10" 

3.125 X 10“ 

6.70 

3.13 

3.04 

0.763 

1.22 

0.49 

H5-WRWH 

HLLR 

HLL 

10® 

3.125 X 10" 

8.04 

3.23 

3.62 

0.747 

0.14 

0.37 

H6-WRWH 

HLLR 

HLL 

10“ 

3.125 X 10" 

5.34 

2.12 

2.31 

0.531 

0.00913 

0.31 

H2-WHWH 

HLL 

HLL 

10" 

3.125 X 10-" 

4.88 

1.47 

1.76 

0.409 

69.6 

0.95 

H3-WHWH 

HLL 

HLL 

10^ 

3.125 X 10“^ 

9.57 

3.54 

4.11 

0.929 

16.1 

1.1 

H4-WHWH 

HLL 

HLL 

10" 

3.125 X 10" 

9.83 

3.93 

4.57 

0.969 

1.77 

0.99 

H5-WHWH 

HLL 

HLL 

10“ 

3.125 X 10" 

9.82 

3.31 

4.57 

0.764 

0.171 

0.95 

H6-WHWH 

HLL 

HLL 

10“ 

3.125 X 10" 

7.41 

2.32 

3.32 

0.573 

0.0125 

0.82 

H2-WRWR 

HLLR 

HLLR 

10" 

3.125 X 10-" 

7.27 

2.15 

2.59 

0.628 

103 

1.9 

H3-WRWR 

HLLR 

HLLR 

10^ 

3.125 X 10-1 

15.8 

5.30 

6.54 

1.40 

25.5 

1.2 

H4-WRWR 

HLLR 

HLLR 

10" 

3.125 X 10" 

24.5 

7.69 

10.4 

1.89 

3.94 

1.9 

H5-WRWR 

HLLR 

HLLR 

10“ 

3.125 X 10" 

27.4 

7.58 

12.0 

1.75 

0.439 

2.3 

H6-WRWR 

HLLR 

HLLR 

10“ 

3.125 X 10" 

133 

24.8 

48.8 

7.31 

0.179 

2.9 
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Table 5: Visco-resistive MHD simulations of the MRI. 


(!) Model 

(2) V X 10^ 

(3) T, X 10" 

{4)/3 

(5) a 

(6) Em X 10" 

(7) Ek X 10" 

(8) wm X 10" 

(9) KJfl X 10" 

(10) a X 10" 

(11) Prrr, 

V3^0-l 

0.5 

0.5 

10" 

3.125 X 10-* 

10.3 

3.76 

4.37 

0.917 

16.9 

1.0 

V4-0-1 

0.5 

0.5 

10* 

3.125 X 10“ 

12.6 

4.18 

5.69 

1.12 

2.18 

1.0 

V5-0-1 

0.5 

0.5 

10" 

3.125 X 10* 

12.5 

3.82 

5.82 

0.979 

0.218 

1.0 

V6^0-l 

0.5 

0.5 

10“ 

3.125 X 10" 

16.1 

4.53 

7.37 

1.21 

0.0274 

1.0 

V3-1-1 

1.0 

1.0 

10^ 

3.125 X 10-* 

8.86 

4.79 

3.65 

0.719 

14.0 

1.0 

V4-1-1 

1.0 

1.0 

10* 

3.125 X 10" 

9.19 

3.04 

4.17 

0.830 

1.60 

1.0 

V5-1-1 

1.0 

1.0 

10" 

3.125 X 10^ 

11.2 

3.13 

5.13 

0.803 

0.190 

1.0 

V6-1-1 

1.0 

1.0 

10“ 

3.125 X 10" 

11.5 

3.26 

5.30 

0.859 

0.0197 

1.0 

V3-2-1 

2.0 

2.0 

10^ 

3.125 X 10-^ 

12.4 

3.70 

4.55 

0.611 

16.5 

1.0 

V4-2-1 

2.0 

2.0 

10* 

3.125 X 10" 

12.3 

2.91 

5.08 

0.822 

1.89 

1.0 

V5-2-1 

2.0 

2.0 

10" 

3.125 X 10* 

16.7 

3.46 

6.79 

0.847 

0.244 

1.0 

V6-2-1 

2.0 

2.0 

10“ 

3.125 X 10" 

16.0 

3.87 

6.48 

0.829 

0.0234 

1.0 

V3-4-1 

4.0 

4.0 

10’ 

3.125 X 10-* 

26.3 

3.03 

8.26 

0.687 

28.6 

1.0 

V4-4-1 

4.0 

4.0 

10* 

3.125 X 10" 

74.2 

7.41 

23.6 

1.72 

8.10 

1.0 

V5-4-1 

4.0 

4.0 

10" 

3.125 X 10* 

216.0 

22.9 

73.2 

6.33 

2.55 

1.0 

V6-4-1 

4.0 

4.0 

10“ 

3.125 X 10" 

358.0 

33.7 

117.0 

9.41 

0.403 

1.0 

V3-2-2 

2.0 

1.0 

10^ 

3.125 X 10-* 

14.9 

3.55 

5.80 

0.907 

21.5 

2.0 

V4-2-2 

2.0 

1.0 

10* 

3.125 X 10" 

12.7 

4.04 

5.49 

0.935 

2.06 

2.0 

V5-2-2 

2.0 

1.0 

10" 

3.125 X 10* 

19.7 

4.25 

8.67 

1.19 

0.316 

2.0 

V6-2-2 

2.0 

1.0 

10“ 

3.125 X 10" 

14.8 

3.12 

6.41 

0.877 

0.0233 

2.0 

V3-4-2 

4.0 

2.0 

10^ 

3.125 X 10-* 

26.4 

3.35 

8.69 

0.924 

30.8 

2.0 

V4-4-2 

4.0 

2.0 

10* 

3.125 X 10" 

48.4 

8.41 

17.3 

2.12 

6.22 

2.0 

V5-4-2 

4.0 

2.0 

10" 

3.125 X 10* 

123.0 

15.4 

41.3 

4.32 

1.46 

2.0 

V6-4-2 

4.0 

2.0 

10“ 

3.125 X 10" 

598.0 

50.8 

188.0 

12.8 

0.642 

2.0 

V3-2-4 

2.0 

0.5 

10^ 

3.125 X 10“* 

21.0 

9.72 

8.38 

1.39 

31.3 

4.0 

V4-2-4 

2.0 

0.5 

10* 

3.125 X 10" 

43.9 

10.4 

18.5 

3.05 

6.90 

4.0 

V5-2-4 

2.0 

0.5 

10" 

3.125 X 10^ 

35.7 

7.07 

15.3 

1.95 

0.551 

4.0 

V6-2-4 

2.0 

0.5 

10“ 

3.125 X 10" 

76.6 

12.8 

30.1 

3.90 

0.109 

4.0 

V3-4-4 

4.0 

1.0 

10’ 

3.125 X 10-* 

28.7 

6.31 

10.1 

1.36 

36.8 

4.0 

V4-4-4 

4.0 

1.0 

10* 

3.125 X 10" 

50.4 

8.20 

19.0 

2.45 

6.86 

4.0 

V5-4-4 

4.0 

1.0 

10" 

3.125 X 10^ 

198.0 

23.9 

65.2 

6.99 

2.31 

4.0 

V6-4-4 

4.0 

1.0 

10“ 

3.125 X 10" 

585.0 

48.0 

183.0 

14.8 

0.633 

4.0 

V3-1-0 

1.0 

2.0 

10" 

3.125 X 10-* 

7.87 

3.34 

3.04 

0.508 

11.4 

0.5 

V4-1-0 

1.0 

2.0 

10* 

3.125 X 10" 

8.62 

2.33 

3.52 

0.573 

1.31 

0.5 

V5M-0 

1.0 

2.0 

10" 

3.125 X 10^ 

7.96 

2.50 

3.44 

0.517 

0.127 

0.5 

V6-1-0 

1.0 

2.0 

10“ 

3.125 X 10" 

8.51 

2.09 

3.49 

0.501 

0.0128 

0.5 







